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14.  ABSTRACT 

This  report  summarizes  the  progress  and  results  for  a  Student  Temporary  Employment  Program  Internship  from  July  2006  through  June 
2007.  The  report  presents  an  introduction  to  the  concept  of  photonic  crystals  through  quantitative  calculations  of  band  gaps  with  the  use  of 
an  open  source  code.  Photonics  is  the  science  and  technology  of  generating,  controlling,  and  detecting  photons,  particularly  in  the  visible 
light  and  near  infrared  spectrum. 

The  first  part  of  the  report  provides  a  detailed  introduction  to  the  principle  of  photonics  from  an  engineering  perspective  and  discusses 
possible  concepts  and  applications  for  the  future.  Sample  calculations  for  phonon  dispersion  curves  are  then  presented  with  the  use  of 
simplified  models  of  crystals  to  demonstrate  universal  principles  of  vibrational  behavior  to  understand  the  properties  of  crystals.  To  this 
end,  a  MATLAB1  program  was  written  to  generate  dispersion  relations  for  two  types  of  reduced  dimensional  lattices. 

The  second  part  of  the  report  demonstrates  the  implementation  of  a  photonic  band  gap  model.  The  calculations  were  perfonned  with  the 
use  of  the  Massachusetts  Institute  of  Technology  photonic  band  gap  code  (MPB).  The  code  was  compiled  and  built  on  the  computers  at 
the  U.S.  Army  Research  Laboratory’s  Major  Shared  Resource  Center,  and  results  verifying  the  installation  are  shown. 

The  report  concludes  with  appendices  detailing  MATLAB  code  for  the  phonon  calculations,  input  decks  for  MPB,  and  an  ancillary 
literature  review  of  ferroelectric  fatigue  in  crystals,  which  was  a  another  topic  of  interest  early  in  the  internship. 
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1.  Background 


Photonics  is  the  study  of  the  creation,  control,  and  detection  of  photons.  The  topic  is  of  increasing 
interest  because  of  projected  value  in  new  types  of  sensor  and  computing  devices.  As  more  is 
learned  about  the  field  and  its  applications,  the  larger  the  support  for  the  research  becomes.  It  is 
closely  related  to  quantum  optics  and  optoelectronics.  Quantum  optics  refer  to  the  principal  research, 
while  photonics  focus  on  the  application  of  the  technology  (7 ).  Opto-electronics  is  the  emerging  field 
of  light-driven  (photon  flow)  electronics  as  opposed  to  the  conventional  electrically  driven  (electron 
flow)  electronics. 

The  basis  of  photonics  is  the  photon1.  The  modern  definition  of  the  photon  was  developed  by 
Albert  Einstein  during  the  early  1900s.  He  originally  called  them  “light  quanta”.  They  were  first 
referred  to  as  photons  by  the  chemist  Gilbert  N.  Lewis  in  1926.  He  based  the  name  on  the  Greek 
word  “phos”  meaning  light  (7).  A  photon  is  the  most  basic  electromagnetic  component.  The 
photon  is  a  quanta  or  “packet”  of  electromagnetic  energy  that  exhibits  particle-like  and  wave-like 
characteristics.  A  photon  is  represented  by  y  or  as  the  product  hv,  h  being  Planck’s  constant  and  v, 
the  frequency  of  the  electromagnetic  wave.  The  momentum  of  a  photon  is  calculated  as  hv/c,  c 
being  the  speed  of  light  which  is  the  constant  velocity  of  all  photons.  Photons  are  massless  and 
have  no  electric  charge.  Photons  are  created  during  several  circumstances  (7): 

1 .  A  charged  particle  undergoes  acceleration; 

2.  An  atomic  particle  drops  from  a  higher  to  lower  energy  state; 

3.  An  atomic  particle  is  destroyed. 

The  field  of  photonics  began  with  the  invention  of  the  laser  in  1960  by  Theodore  Maiman  at 
Hughes  Research  Laboratories.  The  field  began  to  expand  with  the  discovery  of  the  optical  fiber 
in  1970.  The  optical  fiber  functioned  as  a  medium  for  the  transmission  of  light  (7).  The  optical 
fiber  allowed  for  the  controlling  of  photon  flows. 

The  wave-particle  duality  of  light  is  a  crucial  concept  in  the  field  of  photonics.  Einstein  developed 
photons  to  explain  his  experiments  which  showed  that  light  did  not  follow  the  classical  wave 
model.  This  model  depicts  light  as  a  wave  that  propagates  through  a  medium  accordingly.  These 
wave  models  explain  the  refraction,  diffraction,  and  interference  of  light  since  they  are  wave 
behaviors.  On  the  other  hand,  experiments  show  that  photons  do  not  spread  as  they  propagate  and 
they  do  not  split  as  waves  do.  These  are  particle-like  behaviors.  This  comes  to  the  duality  princi¬ 
ple  that  the  photon  particles  produce  the  wave-like  electromagnetic  field.  These  characteristics  of 
light  are  important  factors  in  the  research  and  application  of  photonics. 


1 A  photon  is  a  unit  of  intensity  of  light  at  the  retina  equal  to  the  illumination  received  per  square  millimeter  of  a 
pupillary  area  from  a  surface  having  a  brightness  of  one  candle  per  square  meter. 


1 


2.  Photonic  Crystals 


Photonic  crystals  are  periodic  material  structures  that  affect  the  motion  of  photons  propagating 
through  them.  They  affect  photons  in  the  same  way  a  semiconductor  crystal  affects  the  motion  of 
electrons.  The  crystals  consist  of  periodic  dielectric  structures  that  affect  electromagnetic  wave 
propagation.  A  dielectric  is  a  substance  that  resists  electric  current,  meaning  it  acts  as  an  electrical 
insulator.  In  some  instances,  a  lack  of  substance  or  vacancy  can  act  as  a  dielectric  (e.g.,  air).  The 
photonic  crystals  affect  propagation  by  allowing  electromagnetic  waves  of  certain  wavelengths  to 
pass  while  blocking  others.  A  range  of  blocked  wavelengths  is  called  a  photonic  band  gap  (shown 
in  yellow  in  figure  1).  There  are  some  naturally  occurring  photonic  crystals  such  as  the  gemstone 
opal  and  the  substance  that  comprises  butterfly  wings  (2). 


Figure  1.  Depiction  of  sample  photonic  band  gap.  (Yellow  shaded  region 
shows  a  photonic  band  gap  on  a  photonic  band  diagram.  The 
letters  along  x  axis  represent  the  different  directions  along  the 
crystal  lattice  [5].) 

Diffraction  is  the  basic  principle  behind  the  photonic  crystal  function.  Diffraction  refers  to  dif¬ 
ferent  wave  propagation  characteristics  such  as  bending,  spreading,  and  interference.  Diffraction 
effects  occur  most  when  the  wavelength  and  affecting  medium  structure  are  on  the  same  scale  (2). 
The  visible  spectrum  occurs  on  the  400-  to  700-ntn  wavelength  scale.  As  a  result,  the  feature 
periodicity  of  the  photonic  crystals  must  also  exist  on  this  scale. 

Photonic  crystals  can  affect  propagation  in  one,  two,  and  three  dimensions,  as  shown  in  figure  2. 
The  periodic  structure  and  dielectric  nature  of  the  crystal  detennine  its  ability  to  produce  a  band 
gap.  Lord  Rayleigh  was  the  first  scientist  to  study  the  propagation  of  electromagnetic  waves  in 
periodic  structured  media  in  1887  (2).  He  studied  the  reflective  properties  of  a  crystalline  mineral 
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that  corresponded  to  a  one-dimensional  (1-D)  photonic  crystal.  He  observed  a  small  band  gap 
through  which  light  could  not  propagate  through  the  planes  of  the  crystal. 


Figure  2.  Simplified  representation  of  one-,  two-,  and  three-dimensional  crystals  (4). 


Photonic  crystals  possess  two  types  of  polarizations  by  symmetry:  the  transverse  magnetic  (TM) 
in  which  the  electric  and  magnetic  fields  are  orthogonal  to  one  another,  and  the  transverse  electric 
(TE)  in  which  the  electric  and  magnetic  fields  are  in  the  same  plane.  By  judicious  placement  of 
materials  with  different  high  and  low  indices  inside  an  area  or  volume,  hypothetical  TE  and  TM 
band  gap  materials  may  be  created.  Using  this  simple  principle,  one  can  construct  a  potentially 
infinite  number  of  variations  in  the  topology  of  constituents  to  create  photonic  materials  of  varying 
properties,  as  in  figure  3. 
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Figure  3.  Sample  of  actual  3-D  photonic  band  gap  crystal  structures  ( <D  gives  the 
fdling  fraction  or  percentage  of  unit  cell  occupied  by  material  [5]). 


Eli  Yablonovitch  and  Sajeev  John  both  submitted  independent  papers  on  the  photonic  band  gap  in 
1987,  and  the  search  for  multi-dimensional  band  gap  crystals  began  (6).  Yablonovitch  fabricated 
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the  first  3-D  photonic  band  gap  crystal  while  working  for  Bell  Communications  Research  in  1991 
(6).  The  material  was  composed  of  silicon  and  silicon  dioxide  arranged  in  a  diamond  tetrahedral 
structure.  The  structure  is  now  known  as  yablonovite  (see  figure  4).  Since  their  discovery,  3-D 
band  gaps  have  been  demonstrated  in  simple  cubic  and  face-centered  cubic  structures  as  well. 


Figure  4.  Yablonovite.  (This  was  the  first 

3-D  crystal  with  a  band  gap  that  was 
created  when  holes  were  drilled  into 
a  ceramic  material  [7].) 


3.  Photonic  Crystal  Applications 


The  potential  application  of  photonic  crystals  in  emerging  technologies  is  very  extensive.  The 
volume  of  photonics  research  has  already  produced  many  new  developments  in  the  control  of 
electromagnetic  waves.  One-dimensional  photonic  crystals  are  currently  used  in  thin-film  optic 
applications.  They  are  also  being  used  to  form  reflective  layers  on  mirrors  and  lenses  and  to  create 
color-changing  inks  and  paints.  Two-dimensional  photonic  crystal  fibers  are  produced  by  several 
companies  to  transmit  and  control  light  in  frequencies  that  conventional  fiber  optics  fail  to  transmit 
(6).  Three-dimensional  photonic  crystals  are  much  more  difficult  to  fabricate  and  as  a  result,  have 
not  yet  been  mass  produced  or  applied  to  commercial  products  ( 4 ). 

An  exploratory  simulation  was  done  to  explore  the  use  of  nano  photonic  sensors  for  micro-damage 
detection  (5).  El-Kady  and  Reha  Taha  demonstrated  a  simulation  model  that  is  comprised  of  a  nano 
photonic  sensor  (NPC)  attached  to  a  composite  bar.  When  damage  is  created  in  the  composite  bar, 
the  photonic  sensor’s  band  gap  profile  experiences  a  significant  change.  The  schematic  showing 
this  proposed  setup  is  presented  in  figure  5.  A  slight  change  in  the  NPC’s  surroundings  (damage  in 
the  composite  bar)  will  result  in  a  change  in  the  crystal’s  dimensions.  As  the  crystal’s  dimensions 
change,  so  do  its  optical  properties  (5).  We  can  use  the  NPC  to  quantify  damage  in  the  composite 
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bar  by  monitoring  the  NPC’s  optical  response  before  and  after  damage  is  induced  in  the  bar.  The 
simulation  resulted  in  a  theoretically  verified  damage  metric  (Nano  Spectrum  Index)  that  can  be 
extracted  from  the  frequency  domain  response  of  the  NPC  sensor,  as  seen  in  figure  6  (5).  The 
simulation  will  hopefully  lead  to  the  development  of  NPC  sensors  that  will  allow  for  damage 
detection  at  nano  scales,  which  is  currently  unattainable  with  contemporary  sensors. 


Figure  5.  Schematic  of  simulation  test  setup  (5). 


Figure  6.  Schematic  of  typical  NPC  sensor  response  from  simulation. 

(Frequency  shift  of  NPC  response  is  caused  by  local  damage  in 
the  composite  bar  as  a  result  of  a  non-uniform  strain  [$].) 
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Park  and  Lee  (P)  have  theoretically  investigated  the  effect  of  a  mechanical  stress  applied  to  a 
photonic  crystal.  They  report  a  flexible  photonic  crystal  with  a  tunable  band  structure  that  is  con¬ 
trolled  through  the  use  of  a  nano-/microelectromechanical  system  actuator  (see  figure  7).  For  the 
theoretical  modeling,  the  photonic  crystal  structure  consisted  of  a  pattern  of  high-index  dielectric 
material  and  a  low-index  flexible  polymer  (index  referring  to  the  refractive  index  of  the  material). 
With  the  large  effect  of  physical  change  on  the  band  structure  of  photonic  crystals,  a  large  range  of 
tunability  is  possible  (P).  With  elongations  as  small  as  10%  in  the  photonic  crystal,  a  change  in  the 
refraction  angle  of  an  incident  beam  was  as  large  as  75  degrees  (P).  Figure  8  shows  the  refraction 
angles  for  different  elongations;  these  were  calculated  with  the  theoretically  derived  band  struc¬ 
ture.  The  sizable  change  in  refraction  angle  attributable  to  applied  mechanical  stress  shows 
promise  for  the  use  of  photonic  crystals  as  optical  beam  controlling  devices.  The  goal  of  the 
research  was  to  create  a  system  that  could  control  a  photonic  band  structure  in  real  time.  With 
control  over  the  optical  properties  of  the  photonic  crystal,  the  crystals  can  act  as  optical  switches, 
routers,  and  modulators  on  the  nano-scale  (P). 


Figure  7.  Schematic  of  flexible  photonic  crystal.  (Mechanical  stress  is  applied  to 
the  photonic  crystals  by  micro-electro-mechanical  sensors  actuators  on 
either  side  [9].) 


Figure  8.  Refraction  angle  attributable  to  elongation  (9). 
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Jacobsen  experimentally  found  that  a  linear  electro-optic  effect  can  be  realized  in  silicon  by  break¬ 
ing  the  crystal  symmetry.  The  crystal  symmetry  is  broken  through  the  use  of  a  straining  layer  de¬ 
posited  on  top  of  the  silicon  feature,  e.g.  wave  guide  (see  figure  9).  In  silicon’s  natural  state,  its 
crystal  symmetry  prevents  the  existence  of  a  linear  electro-optic  effect  (10).  As  the  silicon  is 
strained,  its  bulk  refractive  index  fluctuates  in  a  linear  fashion  as  a  function  of  an  externally 
applied  electric  field.  The  experiment  was  conducted  when  a  silicon  nitride  glass  straining  layer 
was  deposited  on  a  silicon  wave  guide;  the  entire  structure  was  on  a  silicon-on-insulator  (SOI) 
wafer  (10).  Application  of  silicon  as  an  opto-electronic  material  is  a  highly  sought-after  objective. 
This  would  allow  for  the  creation  of  combination  electronic  and  optic  components  composed 
entirely  of  silicon.  The  overall  goal  of  the  research  is  to  eventually  use  a  strain-induced  electro¬ 
optic  effect  to  replace  the  electronic  bus  in  computers  with  an  optical  bus  that  performs  at  higher 
speeds  (10). 


Figure  9.  Strain  on  a  crystalline  silicon:  (a)  crystalline  silicon  waveguide  on  an  SOI  wafer; 
(b)  wave  guide  with  a  deposited  straining  layer  (10). 


Many  papers  have  recently  been  published,  documenting  the  advances  in  the  field.  These  advances 
will  hopefully  lead  to  commercialization  of  photonic  crystals  for  advanced  technologies.  Most  of 
the  research  is  in  its  infancy  and  therefore,  many  are  theoretical  simulations  and  analytical  explora¬ 
tions.  As  the  field  and  technology  advances,  more  of  these  applications  will  be  attempted  experi¬ 
mentally.  The  experiments  will  demonstrate  whether  the  applications  are  commercially  viable. 


4.  Dispersion  Curve  MATLAB2  Program 


With  MATLAB,  programs  were  written  to  compute  the  phonon  dispersion  curves  for  different 
materials.  We  adopted  the  use  of  phonons  instead  of  photons  for  didactic  demonstration  of  band 
gaps.  The  code  driving  the  programs  was  written  with  KitteTs  text  (11)  as  a  primary  reference. 

The  text  explained  the  calculation  of  dispersion  curves  for  1-D  crystal  lattices  with  a  monatomic 
structure  as  well  as  with  a  basis.  After  a  code  was  written  to  display  the  dispersion  curves  for  a  1-D 
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crystal,  the  calculations  were  further  expanded  and  a  program  was  created  to  show  the  dispersion 
curves  through  a  2-D  planar  crystal  lattice.  The  MATLAB  program  provides  powerful  graphing 
abilities  that  were  well  suited  to  simulating  the  dispersion  curve  calculations.  The  ensuing  sections 
show  the  derivations  behind  the  calculations,  the  program  code,  and  sample  output  from  the 
programs. 

The  phonon  program  is  a  good  starting  point  for  the  further  study  of  photonics.  The  calculation  of 
the  phonon  dispersion  curves  is  more  pedagogically  simplistic  for  understanding  band  gaps.  The 
phonon  program  provides  a  foundation  for  the  possible  extension  to  graph  photonic  curves.  The 
ability  to  graph  the  frequency  at  which  electromagnetic  energy  quanta  can  pass  through  different 
crystal  lattices  is  an  invaluable  tool.  The  graphs  show  the  band  gaps  and  they  would  be  extremely 
useful  for  the  manufacture  and  control  of  materials  for  photonics  applications. 

The  following  sections  show  the  derivations  of  equations  showing  the  frequency  of  vibration  based 
on  a  wave  vector  k.  These  modes  of  vibration  are  known  as  phonons  and  they  occur  in  structures 
with  rigid  crystal  lattices  (11).  These  phonon  mode  graphs  are  known  as  dispersion  curves.  The 
physical,  such  as  thennal  and  conductive,  properties  of  materials  depend  on  these  phonons. 

The  first  derivation  in  section  4. 1  is  for  the  nonnal  modes  of  a  1-D  monatomic  Bravais  lattice.  In  a 
Bravais  lattice,  all  the  lattice  points  are  the  same.  This  by  default  means  that  the  ions  in  the  crystal 
are  all  the  same  kind,  thus  the  monatomic  label. 

The  derivation  in  section  4.2  is  for  the  normal  modes  of  a  1-D  lattice  with  a  basis.  This  is  still  a 
Bravais  lattice  because  the  basis  is  an  identical  ion. 

Section  4.3  contains  the  derivation  for  the  nonnal  modes  of  a  2-D  monatomic  lattice.  Instead  of 
the  1-D  chain,  this  derivation  solves  for  the  phonon  modes  of  a  2-D  planar  crystal  structure. 

The  next  step  would  be  to  further  develop  the  2-D  code  to  handle  general  crystal  structures.  This 
code  would  need  to  handle  different  crystal  structures  with  different  ions  at  lattice  points  or  struc¬ 
tures  with  ions  that  have  a  basis.  After  that  is  completed,  a  comprehensive  code  for  a  3-D  lattice 
could  be  formulated.  This  would  increase  the  code’s  complexity  immensely  and  would  yield  a 
powerful  computational  tool. 

4.1  Normal  Modes  of  a  One-Dimensional  Monatomic  Bravais  Lattice 

Take  a  line  of  ions  all  of  the  same  mass  M  spaced  at  an  interval  distance  of  a,  as  indicated  in 
figure  10.  This  gives  the  1-D  Bravais  lattice  vectors  as  R  =  na,  where  n  is  a  counting  integral. 
Assuming  the  line  is  envisioned  in  the  horizontal  direction,  the  displacement  of  a  given  ion  from 
its  equilibrium  position  is  u(na).  The  most  basic  way  to  view  the  chain  of  ions  is  as  a  spring-mass 
chain  as  seen  in  figure  10.  This  allows  for  the  spring  potential  energy  equation  to  represent  the 
fundamental  interaction  between  the  ions. 
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Figure  10.  Representation  of  a  1-D  chain  (72). 

U  =  -kx2. 

2 

The  potential  energy  U  is  equated  to  a  function  of  the  spring  constant  k  and  the  ion  displacement. 

k  =  K; 

x  =  \u(na)  -  u([n  +  Y\a)\. 

The  spring  constant  k  is  replaced  with  K=cp  ”(a),  where  <p(x)  represents  the  interaction  energy  of 
two  ions  at  distance  x  apart  along  the  line  (spring  constant).  This  interaction  energy  is  found  with 
the  Lennard- Jones  potential  described  in  appendix  A.  We  find  the  displacement  x  by  taking  the 
displacement  of  the  ion  from  equilibrium  and  subtracting  the  displacement  of  the  next  ion  from  its 
equilibrium.  These  substitutions  give  the  harmonic  potential  energy  equation  as  the  summation 
over  all  ions  in  the  chain: 

jjhann  =  I  ^  [M(Wa)  -  U([n  +  l]a)f  ,  (1) 

2  n 

Qjjharm 

Mii(na )  = - =  -K[2u(na)  -  u([n  -  l]n)  -  u([n  +  l]a)] .  (2) 

du(na ) 

The  equations  of  motion  of  the  ions  show  the  theoretical  interaction  of  two  neighboring  ions 
connected  by  a  perfect  spring  (11).  We  find  the  motion  equation  by  taking  the  derivative  of  the 
harmonic  potential  energy  with  respect  to  the  ion’s  displacement. 

If  the  chain  is  considered  a  finite  system  with  N  ions,  boundary  conditions  can  be  implemented  to 
retrieve  solutions  from  the  motion  equation.  In  order  to  apply  boundary  conditions,  the  finite 
chain  must  be  viewed  as  a  loop  so  the  N+l  ion  is  the  first  ion  of  the  chain  (11).  These  boundary 
conditions  are  expressed  mathematically  as 

u([N  +  l]n)  =  u(a);  u( 0)  =  u(Na) .  (3) 

Solutions  to  the  equation  of  motion  in  equation  2  take  the  form 

u(na,t)  oc  ei(kna-M) .  (4) 

Applying  the  boundary  conditions  of  equation  3  leads  to  the  wave  vector  k  equation  to  be 
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(5) 


ikNa  |  j  2jt  H 

e  =1;  k  = - . 

a  N 

Only  N  values  of  k  consistent  with  equation  5  yield  distinct  solutions,  taken  between  -n/a  and  n/a 
in  order  to  view  the  dispersion  curve  about  the  zero  wave  vector.  This  range  covers  2n/a\  any  of 
the  distinct  solutions  shifted  by  2n/a  will  return  the  same  result,  meaning  that  this  range  produces 
the  entire  distinct  dispersion  curve. 

Substituting  the  solutions  of  the  form  in  equation  4  into  the  equation  2  of  motion  gives  the 
following: 


-Mco~e 


2  i(kna-cot)  _ 


ika  ika  ~i  i(kna-cot) 


-K\2  -  e  -e'Ka]e 


Euler’s  formulas  (subtraction)  can  be  used  to  simplify  equation  6  to  a  more  manageable  and 
graphing  friendly  result.  The  Euler  subtraction  is  shown  below: 

e'x  =  cos  x  +  i  sin  x 

-  =  cos  x  -  i  sin  x  =>  2  cos  x  =  e~n  +  ea 

cos  x  +  cos  x  -  =  e,x 

-K[ 2  -  eika  -  eika  ]  =>  — 2AT[1  - 1  cos  ka] . 


Using  the  simplification  in  equation  6, 


-M(Q2ei(km,-6,t)  =  -2K(\  -  cos  ka)el(kna~<ot) , 
-Mor  =  -2K(\- cos  ka), 

U,  f2^0  -  cos  ka) 


co(k)  = 


The  final  equation  of  the  dispersion  curve  frequency  as  a  function  of  the  wave  vector  k  is  shown 
by  equation  8;  it  is  found  with  the  half  angle  to  further  simplify  result  in  equation  7: 

.  2  l-cos(2x)  _  .  2  r  1  ^  , 

sin“x  = - =>2sin  — x  =l-cosx, 

2  v  2  J 


co(k)  = , 


2K(l  -cos  ka) 


2 K  2  sin  —  ka 
{  2 


„  K  .1, 
co(k)  =  2,  —  sin  —  ka  . 
VM  2 


Solutions  describing  the  actual  displacement  of  the  ions  are  given  by  the  real  or  imaginary  parts 
of  equation  4,  i.e., 
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u(na,t )  oc 


cos(kna  -  cot) 
sin (kna  -  cot) 


(9) 


4.2  Normal  Modes  of  a  One-Dimensional  Lattice  with  a  Basis 

The  normal  mode  calculation  of  a  1-D  lattice  with  a  basis  is  very  similar  to  the  monatomic  proce¬ 
dure.  The  representation  in  figure  1 1  shows  the  physical  schematic  as  well  as  the  spring-mass 
chain  visualization.  The  ions  all  have  the  same  mass,  M.  In  this  crystal  structure,  there  are  two 
equilibrium  positions,  na  and  na+d,  for  the  ion  and  its  basis,  respectively.  It  is  assumed  that  only 
the  closest  neighbors  interact. 

—"I  a—d  — - 1  d  H—  — H  a  h— 

•  X  •  X*  X*  X*  X  •  X*  X*  X* 

na  (n  +  \)a  (n  +  2)a  (n  +  3)a  (n  +  4)a  (n  +  5)a  (n  +  6)a  (n+l)a 

'WMP  G-spring 

AAA,  A-spring _ 

Figure  11.  Representation  of  a  1-D  chain  with  a  basis  ( 11 ). 


The  spring-mass  chain  representation  can  be  expressed  mathematically  as 


U  =  —Lx ,2  +  —k2x12 . 
2  2 


The  potential  energy  is  equal  to  the  spring  constant  and  the  ion  displacement.  For  the  lattice  with 
a  basis,  both  ions  must  be  accounted  for  and  because  of  the  difference  in  separation  distances,  the 
different  spring  constants  and  displacements  must  be  recognized. 

kx=K\  k2  =  G, 

Xj  =  i/j(Ha)-u2(no), 

x2  =u2(na)-ul([n  +  l]a)\, 

-ui(na)  is  the  displacement  of  ion  oscillating  around  site  na, 

-U2(na)  is  the  displacement  of  ion  oscillating  around  site  na+d. 

As  before,  the  spring  constant  k  is  replaced  with  K,  G=cp  ”(a),  where  <p(x)  represents  the  interaction 
energy  of  two  ions  at  distance  x  apart  along  the  line  (spring  constant).  This  interaction  energy  is 
found  with  the  Lennard-Jones  potential  (appendix  A).  These  substitutions  give  the  hannonic 
potential  energy  equation  as  the  summation  over  all  ions  in  the  chain: 

jjhann  =  (77a)  - u2 (na)]2  +  —  ^ [u2 (na)  - u1  (\n  +  l]a)]2  . 

^  n  ^  n 
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From  the  harmonic  potential  energy  equation,  the  equations  of  motion  can  be  found  if  we  take 
the  second  derivative: 

^j^jharm 

Mil,  ( na )  = - 

dux(na) 

=  -K [ux ( na )  - u2 (na)]  -  G[ux (na)  - u2([n  - l]a)], 

^jj-harm 

Mii ,  (na)  = - 

du2(na) 

=  -K \u2  (na)  -  ux  (na)]  -  G[ii2  (na)  -  ux([n  +  l]n )] . 

There  will  be  two  solutions  to  the  motion  system  of  equations.  This  is  because  there  will  be  two 
modes  for  the  ion  and  its  basis.  They  can  be  oscillating  at  different  frequencies  about  their 
respective  equilibrium  positions.  These  solutions  will  take  the  form 

ux(na)  =  £xel(kna-mt) , 
u2  (na)  =  £2eI(kna~mt). 

The  values  of  Si  and  82  are  ratios  that  will  specify  the  relative  amplitude  and  phase  of  vibration  of 
the  ions  within  each  lattice  cell  (11).  The  same  boundary  conditions  used  for  the  1-D  monatomic 
chain  are  used  again  to  retrieve  solutions  from  the  equations  of  motion.  Substituting  the  solution 
equations  into  the  equations  of  motion  produces  two  coupled  equations  by  canceling  the  common 
exponential  factor,  leading  to 

[Mffl2  ~(K  +  G)]£x  +(K  +  Ge~ika  )e2  =  0,  and 
(K  +  Geika)£x  +  [Mar  ~(K  +  G)]e2  =  0. 

The  homogeneous  equations  will  have  a  solution,  provided  the  determinant  of  the  coefficients 
vanishes,  leaving 

[Mco2  -( K  +  G )]2  =\K  +  Ge~ika\2  =  K2  +  G2  +  2KG cos ka  . 

This  equation  is  satisfied  by  two  positive  values  of  ro.  These  two  values  correspond  to  two 
different  modes  on  the  dispersion  diagram  and  are  found  from  the  roots  of 

or  =  ^^-±  —  ^K2  +G2  +2KG  cos  ka  . 

M  M 

4.3  One-Dimensional  MATLAB  Results 

The  MATLAB  code  can  be  modified  to  find  the  1-D  phonon  normal  modes  for  a  monatomic  chain 
as  well  as  a  chain  with  a  basis.  The  code  in  appendix  B  solves  for  the  1-D  normal  modes.  In  that 
code,  the  two  solutions  to  0  (ion  oscillation  frequency)  can  be  seen  in  the  wave  vector  loop. 
Including  both  of  those  equations  gives  the  solution  for  a  1-D  chain  with  a  basis  (see  figure  12). 

If  one  of  the  equations  is  removed,  the  dispersion  curve  shows  the  normal  mode  for  a  1-D  chain 
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without  a  basis  (see  figure  13).  The  dispersion  relations  are  symmetric  about  the  zero  wave  vector. 
As  mentioned  in  the  derivation  for  the  code  in  section  4.1,  the  wave  vectors  from  -n/a  to  n/a  gene¬ 
rate  the  entire  curve  for  the  dispersion  relation.  The  curve  would  simply  be  repeated  if  a  larger 
range  of  k  were  graphed.  At  small  values  of  the  wave  vector,  the  relation  is  linear;  the  wave  vector 
is  short  enough  to  be  comparable  to  the  particle  spacing;  thus,  the  linearity  ceases  and  the  disper¬ 
sion  curve  flattens  as  the  group  velocity  tends  to  zero  (11).  The  group  velocity  is  the  velocity  at 
which  the  change  in  the  shape  of  the  wave's  amplitude  propagates  through  a  medium  (1).  This 
corresponds  to  the  slope  of  the  dispersion  curve. 


3.5 
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i - i - i - i  i  i  — 

- Optical  Branch 

- Acoustic  Branch 


2.5 
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0.5 


O' - <- 
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0.2  0.4 
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Figure  12.  Dispersion  curve  for  a  1-D  chain  with  a  basis. 


Figure  13.  Dispersion  curve  for  a  monatomic  1-D  chain. 
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The  two  functions  of  frequency,  to,  are  known  as  the  two  branches  of  the  dispersion  relation  (II). 
The  lower  branch  is  known  as  the  acoustic  branch  because  its  dispersion  curve  is  in  a  form  similar 
to  sound  waves.  The  upper  branch  is  referred  to  as  the  optical  branch  because  the  longer  wave¬ 
length  optical  modes  in  crystals  can  interact  with  electromagnetic  radiation.  The  longer  wave¬ 
length  modes  are  also  responsible  for  the  optical  characteristics  of  the  crystal  (11). 

Figure  14  shows  the  vibration  of  the  ions  in  the  1-D  chain  with  a  basis  for  the  acoustic  and  optical 
branch.  If  the  wave  vector  is  taken  to  be  zero,  the  frequency  is  equal  to  zero.  This  causes  the  ion 
and  basis  to  have  the  same  amplitude  and  phase.  The  pair  moves  in  phase  as  the  acoustic  mode 
(2).  For  the  optical  mode,  the  ion  and  its  basis  move  at  the  same  amplitude  but  out  of  phase;  this 
results  in  the  pair’s  center  of  mass  remaining  constant.  This  optical  mode  occurs  with  a  frequency 
of  vibration  in  the  infrared  region,  thus  the  optical  branch  label. 


4.4  Normal  Modes  of  a  Two-Dimensional  Monatomic  Lattice 


To  calculate  the  normal  modes  of  a  2-D  lattice,  the  methods  for  the  1-D  chain  were  modified  and 
expanded.  We  accomplished  this  by  taking  the  planar  lattice  as  a  mesh  of  three  chains:  a  hori¬ 
zontal,  a  vertical,  and  a  diagonal.  As  shown  in  figure  15,  the  lattice  has  an  arbitrary  spacing  of  a, 
b,  and  the  diagonal  distance  in  which  those  result.  For  simplicity,  the  lattice  is  made  to  be  mona¬ 
tomic.  The  1  and  n  variables  serve  as  counters  to  designate  which  ion  in  the  lattice  is  being  used. 


Applying  these  three  chains  results  in  an  equation  for  the  potential  energy,  given  by 


u  =  1k^[u 
£  1 


J  +  joy  [»,„  . 

"  n 


Nearest  Neighbors 


l  n 

Next  Nearest  Neighbors 
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Figure  15.  Representation  of  a  2-D  monatomic 
lattice. 


As  before,  we  find  the  equation  of  motion  for  the  transverse  modes  through  the  medium  by  taking 
the  second  derivative  of  the  potential  energy  equation.  K,  G,  and  H  are  the  three  “spring  constants” 
or  interaction  energy  between  the  ions.  These  values  were  found  with  the  Lennard-Jones  potential. 
For  the  2-D  modes,  we  made  the  calculation  more  accurate  by  also  looking  at  the  next  closest 
neighbors  for  the  interaction  energy.  If  the  lattice  spacing  constants  a  and  b  are  equal  (square 
lattice),  the  values  of  K  and  G  will  be  the  same  because  of  the  monatomic  stipulation.  The  resulting 
equation  of  motion  is 

MUi,n  =  ~K  [2w/,„  -  u,-x.n  -  UM ]  -  G  [2w/,„  -  1  -  u,,,+ 1  ]  -  2H  [2u,,n  -  K/-U-1  -  W/+1,„+1  ] 

With  a  monatomic  lattice,  there  is  only  one  solution  form  to  satisfy  the  equation  of  motion,  which 
is 

i(krla+k„nb-(Ot) 

UUn  =e 

Substituting  the  solution  into  the  equations  of  motion  and  solving  for  the  frequency,  co  is  shown 
below: 


arM  =  -K  [2  -  eKa  -  eik*a  ]  -  G 


m  ~ikvb  ikvb 

2  —  e  —  e  - 


-2  H 


^  - i(kYa+kvb )  i(kra+kvb ) 

2  —  e  —e 


arM  =  -2 K  [l  -  cos(kYa)]  -  2G  [l  -  cos(£vh)]  -  4 H  [l  -  cos (kxa  +  kyb )] 


ao  =  J— 


M 


[if  (1  -  cos (kxa))  +  G(1  -  cos (kyb))  +  2if  (1  -  cos {kxa  +  kyb))~j . 


The  resulting  normal  mode  equation  is 


9  \K 

co  =  2J — 
\M 


■  Xi 

sin —ka 
2  x 


+  2.1— 

M 


sin~k  b 

2  -v 


.  ,H 
+  4J — 
M 


sin  (~^kxa  +  kyb) 
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With  the  2-D  planar  lattice,  it  is  possible  to  have  wave  vectors  in  multiple  directions  as  opposed  to 
the  1-D  chain  which  could  only  support  a  wave  vector  in  the  chain’s  direction.  These  allowable 
wave  vectors  are 


r  2n  l  „  In  n  „ 

k  = - xH - v: 

a  L  b  N 


l,n  =  0,±  1,±  2...  + 


L,  N 

~Y~ 


4.5  Two-Dimensional  MATLAB  Results 

The  MATLAB  code  for  the  2-D  planar  lattice,  as  shown  in  appendix  C,  was  more  complicated  than 
the  1-D  chain.  The  main  complexity  came  from  the  wave  vector  consideration.  In  a  1-D  chain,  the 
only  possible  wave  vector  is  the  vector  traveling  through  the  chain.  In  the  2-D  plane,  there  may 
exist  horizontal,  vertical,  and  diagonal  wave  vectors.  Figure  16  depicts  the  dispersion  relation  for  a 
wave  vector  with  only  a  horizontal  component  (the  vertical  wave  vector  component  was  set  to  zero). 

kx=0 


Figure  16.  Dispersion  curve  for  a  2-D  lattice.  (The  wave  vector  is  only 
in  the  horizontal  direction.) 


In  figure  17,  the  horizontal  and  vertical  wave  vector  components  were  set  to  equal  each  other  so 
the  wave  travels  diagonally  through  the  planar  lattice.  This  results  in  a  3-D  graph  since  the 
dispersion  relation  has  components  in  the  vertical  and  horizontal. 

kx—ky 

In  figures  16  and  17,  it  is  clear  that  the  depicted  dispersion  relation  only  has  an  acoustic  mode. 
These  curves  are  very  similar  to  what  was  seen  with  the  1-D  monatomic  chain.  With  the  planar 
lattice  consisting  of  only  a  single  type  of  ion,  it  is  expected  there  would  be  no  optical  mode.  If 
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another  ion  were  introduced  into  the  lattice  or  a  basis  were  added  to  each  lattice  point,  an  optical 
branch  would  appear  in  the  dispersion  relation. 


Figure  17.  Dispersion  curve  for  a  2-D  lattice.  (The  wave  vector  is  diagonal.) 


5.  Massachusetts  Institute  of  Technology  (MIT)  Photonic  Bands 


The  MIT  photonic  bands  (MPB)  package  is  a  Unix3-based  program  for  computing  band  structures 
of  periodic  structures.  The  package  was  developed  by  Steven  G.  Johnson  and  the  Joannopoulos 
Ab  Initio  Physics  group  at  MIT.  The  program  computes  the  eigenstates  (harmonic  modes)  of 
Maxwell’s  equations  in  periodic  structures,  given  arbitrary  wave  vectors  (13).  The  program  uses 
vector-based  math  and  3-D  analytical  methods. 

The  program  is  a  large  step  above  the  MATLAB  code  introduced  in  the  previous  section.  It  was 
a  logical  step  in  the  research  progression;  the  MPB  package  works  well  with  photonic  crystals 
whereas  the  MATLAB  code  was  only  applicable  to  phonon  mode  calculations.  The  MPB  package 
can  accept  a  wide  range  of  periodic  structures  and  calculates  the  corresponding  dispersion  relations. 

The  MPB  package  runs  a  set  of  provided  specifications  that  are  set  in  a  control  file.  The  control 
files  are  formatted  text  documents  of  the  fonn  “*.ctl”.  These  control  files  set  the  geometry  of  the 
structure,  the  number  of  eigenvectors  to  compute,  what  values  and  information  to  output,  and  any 
other  important  specifications.  When  run,  the  control  file  activates  and  accesses  a  library  (libctl) 
of  programming  commands  that  tell  the  MPB  program  what  to  compute.  This  allows  the  creation 
of  the  control  files  to  be  very  simple.  To  run  the  MPB  program,  a  command  such  as  the  following 
is  typed  in  the  Unix  prompt: 


Unix  is  a  registered  trademark  of  The  Open  Group. 
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unix%  mpb  *.ctl  >&  *.out 


This  command  tells  the  MPB  program  to  run  the  *.ctl  file  and  output  the  results  in  the  *.out  file. 
This  *.out  file  is  a  readable  text  file  that  will  provide  the  results  of  the  calculations  the  program  has 
run.  Several  commands  that  can  be  invoked  at  the  prompt  provide  specific  information  contained 
in  the  *.out  without  the  user  opening  the  whole  document. 

•  unix%  grep  Gap  *.out 

This  command  will  provide  any  band  gaps  in  the  computed  *.out  file.  It  is  important  to  note  that 
any  gaps  <1%  are  most  likely  attributable  to  band  crossings  and  are  not  true  band  gaps;  they  have 
been  shown  in  red  in  the  sample  output  below. 

Gap  from  band  1  (0.275065617068082)  to  band  2  (0.446289918847647),  47.4729292989213% 

Gap  from  band  3  (0.563582903703468)  to  band  4  (0.59305906621551 1),  5.0968516236891% 

Gap  from  band  4  (0.791 161222813268)  to  band  5  (0.79204273 1370125),  0.1 11357548663006% 

Gap  from  band  5  (0.838730315053238)  to  band  6  (0.840305955160638),  0.187683867865441% 

Gap  from  band  6  (0.869285340346465)  to  band  7  (0.873496724070656),  0.483294361375001% 

Gap  from  band  4  (0.821658212109559)  to  band  5  (0.864454087942874),  5.07627823271133% 


•  unix%  grep  tmfreqs  *.out  >  *.tm.dat 
unix%  grep  tefreqs  *.out  >  *.te.dat 

This  pair  of  commands  takes  the  TE  and  TM  bands  from  the  *.out  output  file  and  organizes  the 
band  data  into  a  unifonn  list.  This  makes  it  very  easy  to  take  the  band  structure  data  and  use  it  for 
whatever  is  desired.  A  sample  *.te.dat  file  is  shown  below. 

tefreqs:,  1,  0,  0,  0,  0,  0,  1,  1,  1,  1,  1.41421,  1.41421,  1.41421 

tefreqs:,  2,  0.1,  0,  0,  0.1,  0.1,  0.9,  1.00499,  1.00499,  1.1,  1.34536,  1.34536,  1.48661 

tefreqs:,  3,  0.2,  0,  0,  0.2,  0.2,  0.8,  1.0198,  1.0198,  1.2,  1.28062,  1.28062,  1.56205 

tefreqs:,  4,  0.3,  0,  0,  0.3,  0.3,  0.7,  1.04403,  1.04403,  1.22066,  1.22066,  1.3,  1.64012 

tefreqs:,  5,  0.4,  0,  0,  0.4,  0.4,  0.6,  1.07703,  1.07703,  1.16619,  1.16619,  1.4,  1.72047 

tefreqs:,  6,  0.5,  0,  0,  0.5,  0.5,  0.5,  1.11803,  1.11803,  1.11803,  1.11803,  1.5,  1.80278 

tefreqs:,  7,  0.5,  0.1,  0,  0.509902,  0.509902,  0.509902,  1.02956,  1.02956,  1.2083,  1.2083,  1.50333,  1.74929 

tefreqs:,  8,  0.5,  0.2,  0,  0.538516,  0.538516,  0.538518,  0.943398,  0.943399,  1.3,  1.3,  1.51327,  1.7 

tefreqs:,  9,  0.5,  0.3,  0,  0.583095,  0.583095,  0.583095,  0.860233,  0.860233,  1.39284,  1.39284,  1.52971,  1.52971 

tefreqs:,  10,  0.5,  0.4,  0,  0.640312,  0.640312,  0.640312,  0.781025,  0.781025,  1.48661,  1.48661,  1.55242,1.55242 

tefreqs:,  11,  0.5,  0.5,  0,  0.707107,  0.707107,  0.707107,  0.707107,  0.707107,  1.58114,  1.58114,  1.58114,1.58114 

tefreqs:,  12,  0.4,  0.4,  0,  0.565685,  0.565685,  0.72111,  0.72111,  0.848528,  1.45602,  1.45602,  1.52315,  1.64924 

tefreqs:,  13,  0.3,  0.3,  0,  0.424264,  0.424264,  0.762898,  0.763656,  0.990966,  1.33459,  1.33471,  1.47686,1.72306 

tefreqs:,  14,  0.2,  0.2,  0,  0.282843,  0.282843,  0.824621,  0.824621,  1.13137,  1.21655,  1.21655,  1.44222,  1.44222 

tefreqs:,  15,  0.1,  0.1,  0,  0.141421,  0.141421,  0.905539,  0.905539,  1.10454,  1.10454,  1.27279,  1.42127,  1.42127 

tefreqs:,  16,  0,  0,  0,  0,  0,  1,  1,  1,  1,  1.41421,  1.41421,  1.41421 

The  following  three  sections  demonstrate  just  a  small  sample  of  the  capabilities  of  the  MPB 
program.  Sections  5.1  and  5.2  show  the  analyses  of  a  diamond  and  tri-rod  structure,  respectively. 
Section  5.3  shows  how  a  new  structure  can  be  defined  and  analyzed  with  the  control  file. 
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To  construct  graphs  from  the  output  files  (*.out),  the  following  procedure  was  followed.  It  is 
important  to  note  that  all  the  steps  can  be  done  with  Unix-based  programs,  but  for  ease  of  use,  the 
data  were  brought  into  the  Windows4  OS  for  manipulation. 

1 .  Run  the  following  to  create  an  output  file  containing  the  data  for  the  desired  periodic  structure: 
unix%  mpb  *.ctl  >&  *.out 

2.  Run  the  following  to  organize  the  output  file  data  into  a  manageable  format: 
unix%  grep  tmfreqs  *.out  >  *.tm.dat 

unix%  grep  tefreqs  *.out  >  *.te.dat 

3.  With  Excel5  the  data  can  be  imported  from  the  *.dat  file  into  the  spreadsheet.  It  is  crucial 
in  this  step  to  use  comma  delimitation  to  have  to  data  organized  in  the  spreadsheet  to 
greatly  accelerate  the  graphing  process. 

4.  In  Excel,  the  data  can  be  easily  graphed  and  analyzed. 

5,1  Diamond  Structure 

The  MPB  package  came  with  a  folder  of  sample  structures.  This  folder  contained  the  control  files 
for  a  range  of  different  periodic  structures  to  demonstrate  the  MPB  software’s  capabilities.  The 
on-line  MPB  tutorial  (13)  provides  a  walk-through  of  the  diamond  structure  control  file;  the  walk¬ 
through  is  covered  in  this  section.  The  diamond  control  file  (diamond. ctl)  sets  up  the  periodic 
structure  of  a  3-D  diamond  lattice  of  dielectric  spheres  in  air.  A  3-D  model  of  this  structure  is 
shown  in  figure  18.  The  control  file  is  given  in  appendix  D.  This  structure  and  corresponding  full 
band  gap  were  first  introduced  by  Ho,  Chan,  and  Soukoulis  (14). 

We  calculated  the  bands  by  running  the  following  command: 
unix%  mpb  diamond. ctl  >&  diamond. out 

Based  on  Ho,  Chan,  and  Soukoulis  (14),  a  full  band  gap  should  exist  in  the  band  structure.  This  is 
explored  by  the  following  command: 

unix%  grep  Gap  diamond. out 

Gap  from  band  2  (0.396348703007373)  to  band  3  (0.440813418580596),  10.6227251392791% 

This  output  shows  that  there  is  a  10.6%  band  gap  from  band  2  to  band  3.  This  gap  can  be  seen  in 
figure  19  from  the  tutorial  and  figure  20  which  was  generated  with  Excel  to  graph  the  MPB  output 
data.  Figures  18  and  19  are  identical,  which  verifies  that  the  MPB  code  and  resulting  data  files 
were  used  correctly.  Figure  19  also  matches  a  figure  from  the  paper  (14).  The  figure  from  the 
paper  was  generated  by  other  means;  it  is  encouraging  that  it  matches  the  MPB  output. 


^Windows  is  a  trademark  of  Microsoft  Corporation. 
5Excel  is  a  trademark  of  Microsoft  Corporation. 
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Figure  18.  Diamond  structure  depiction  (13)  (MIT 
photonic-band  tutorial  figure). 


Figure  19.  Diamond  band  diagram  (14)  (MIT  photonic-band  tutorial  figure). 
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Figure  20.  Diamond  band  diagram  (MIT  photonic-band  program  generated  plot). 


5.2  Tri-rods  Structure 

The  second  periodic  structure  MPB  example  that  was  covered  in  the  on-line  MPB  tutorial  (13)  was 
the  tri-rods  structure  control  file  (tri-rods. ctl).  The  walk-through  from  the  on-line  tutorial  is  covered 
in  this  section.  The  structure  is  a  2-D  triangular  lattice  of  dielectric  rods  (figure  19).  The  control 
file,  tri-rods. ctl  is  shown  in  appendix  E. 


Figure  21.  Rods  in  a  triangular 
lattice  (75),  top  view. 

We  calculated  the  bands  by  running  the  following  command: 

unix%  rnpb  tri-rods. ctl  >&  tri-rods. out 
The  band  gaps  were  verified  with  the  following  command: 
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unix%  grep  Gap  tri-rods. out 

Gap  from  band  1  (0.275065617068082)  to  band  2  (0.446289918847647),  47.4729292989213% 

Gap  from  band  3  (0.563582903703468)  to  band  4  (0.593059066215511),  5.0968516236891% 

Gap  from  band  4  (0.791161222813268)  to  band  5  (0.792042731370125),  0.111357548663006% 

Gap  from  band  5  (0.838730315053238)  to  band  6  (0.840305955160638),  0.187683867865441% 

Gap  from  band  6  (0.869285340346465)  to  band  7  (0.873496724070656),  0.483294361375001% 

Gap  from  band  4  (0.821658212109559)  to  band  5  (0.864454087942874),  5.07627823271133% 

The  gaps  in  red  are  less  than  1%  and  are  therefore  bands  crossing  one  another.  The  other  three 
bands  are  significant  in  size  and  can  be  seen  on  the  band  structure  diagram  (see  figure  22).  It  is 
important  to  note  that  unlike  the  diamond  structure,  there  is  no  complete  band  gap.  There  are  only 
partial  TE  and  TM  band  gaps.  In  figure  22,  the  light  blue  area  represents  TM  band  gaps  and  the 
light  red  represents  the  TE  band  gap. 

Figure  23  shows  the  Excel  graph  generated  from  the  output  ( tri-rods.  *.dat).  It  can  be  seen  that  the 
band  structure  in  figure  23  identically  matches  the  band  structure  diagram  in  figure  22  from  the  on¬ 
line  MPB  tutorial. 


Figure  22.  Tri-rods  band  diagram  (13)  (MIT  photonic  band  tutorial  figure). 
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Figure  23.  Tri-rods  band  diagram  (MIT  photonic  band  program-generated  plot). 


5.3  Sample  Structure 

After  the  two  examples  are  complete,  the  on-line  MPB  tutorial  (13)  shows  how  to  create  a  control 
file.  This  is  shown  in  the  user  tutorial  section.  After  one  becomes  familiar  with  the  basic  commands 
through  the  user  tutorial,  the  user  reference  section  on  line  (13)  can  be  used  to  further  explore  com¬ 
mands  that  can  be  used  in  the  control  file.  To  test  the  program  and  control  file  function,  a  sample 
structure  was  developed  and  run.  The  sample  control  file  (sample. ctl)  is  shown  in  figure  24.  The 
control  file  is  very  short  and  concise.  This  is  part  of  the  advantage  of  the  library  (libctl)  discussed 
before.  It  allows  short  commands  such  as  "set!  *"  to  call  up  a  longer  code  in  the  library  that  the 
MPB  program  then  executes.  The  sample  structure  is  a  3-D  cube  of  dielectric  material  with  an  air 
hole  at  its  center. 

We  calculated  the  bands  for  the  sample  structure  by  running  the  following  command: 
unix%  mpb  sample. ctl  >&  sample. out 

To  look  at  the  band  gap,  the  familiar  “unix%  grep  Gap  sample. out”  was  used,  but  no  output  was 
provided.  This  is  because  there  are  no  partial  or  full  band  gaps  in  the  structure.  This  is  evident 
in  the  band  diagram  in  figure  25. 
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(set!  num-bands  8) 

(set!  k-points  (list  (vector3  0  0  0)  ;  Gamma 

(vector3  0.5  0  0)  ;  X 
(vector3  0.5  0.5  0) ;  M 
(vector3  0  0  0)))  ;  Gamma 

(set!  k-points  (interpolate  4  k-points)) 

;  A  unit  cube  of  dielectric  material  with  a  spherical  air  hole  of  radius  0.2  at 
;  its  center, 

(set!  geometry  (list 

(make  block  (center  12  3)  (material  (make  dielectric  (epsilon  12)))  (size  111)) 
(make  sphere  (center  12  3)  (material  air)  (radius  0.2)))) 

(set!  geometry-lattice  (make  lattice  (size  1  1  no-size))) 

(set!  resolution  32) 

(run-tm) 

(run-te) 


Figure  24.  Sample  control  file  for  MPB. 


Figure  25.  Sample  band  structure  (MIT  photonic  band  program-generated  plot). 
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5.4  MIT  Photonic-Bands  Conclusion 

The  MPB  package  is  a  powerful  program  to  use  in  the  study  of  photonic  crystals.  After  the  pre¬ 
liminary  study  of  the  background  of  photonics  and  photonic  crystals,  the  program’s  capabilities  are 
easy  to  appreciate.  The  ability  to  specify  the  geometry  of  periodic  structures  with  easy  commands 
in  the  control  file  allows  for  the  fast  study  of  different  structures.  The  MPB  program  is  useful  as  a 
concept  evaluator  to  analyze  and  predict  the  band  structures  of  conceived  designs  through 
computational  modeling  before  we  perform  physical  experiments. 

The  MPB  package  proved  to  be  a  useful  educational  tool.  The  control  file  feature  allows  the 
program  to  be  very  user  friendly.  The  basic  commands  used  in  the  control  file  can  be  quickly 
learned  from  the  on-line  tutorial  and  user  reference,  leading  to  immediate  results  of  photonic 
properties  of  materials. 
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Appendix  A.  Lennard-Jones  Potential 


The  Lennard-Jones  potential  is  a  mathematical  fonnula  that  approximates  the  interaction  energy 
between  two  ions,  based  on  their  radius  of  separation  r.  As  seen  in  figure  A- 1 ,  at  close  range,  there 
is  a  repulsive  energy  that  goes  to  infinity.  Farther  away,  there  is  an  attractive  force  that  tends 
toward  zero  as  the  radius  of  separation  is  increased.  In  the  derivation  below,  &(r)  is  the  value  of 
the  interaction  potential.  The  derivative  of  &(r)  is  taken  twice  to  find  the  “spring  constant”  or 
interaction  energy  between  two  ions.  This  interaction  energy  was  used  in  the  nonnal  phonon  mode 
derivations  and  MATLAB  programs. 

f(r)  =  K,G,H 
r  =  a,b, 

(f j>(r )  =  4 


<j)(r)  =  4 a  a12 


d2(/){r ) 
dr2 


</>"(r )  =  8^ 


78cr12 

r14  J 
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Figure  A- 1.  Lennard- Jones  potential.  (The  graph  shows  how  closely  the 

Lennard-Jones  potential  fits  the  empirical  curve  for  argon  atoms. 
The  x-axis  is  the  radius  of  separation  in  angstroms  [7].) 
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Appendix  B.  MATLAB  Code:  Dispersion  Curves  for  One-Dimensional 
Lattice 


clear  all 
clc 

%function  omega  =  dispersion(phi) 
syms  L  N  omega  k  kx  ky  a  b  K  G  M 

%Test  Input 
L=10; 

N=10; 

a=4.33; 

b=20; 

K=l; 

G=2; 

M=l; 

j=i; 

%Set  loop  to  cover  entire  lattice 
e=-N/2; 
r=N/2; 
for  i  =  e: .  1  :r 
%Wave  vector 
kx(j)=  (2*pi()/a)*(i/L); 

%Dispersion  Curve  Equation 

omega  1  (j)  =  sqrt(((K+G)/M)+(  1  /M)*  sqrt(KA2+GA2+2  *K*  G  *  cos(kx(j)  *  a))); 
omega2(j)  =  sqrt(((K+G)/M)-(l/M)*sqrt(KA2+GA2+2*KH5G*cos(kx(j)*a))); 

j=j+i; 

end 

kx=double(kx); 

omega  1  =double(omega  1 ); 

omega2=double(omega2); 

plot(kx,  omega  1,  kx,  omega2) 

legend(’Optical  Branch'/Acoustic  Branch','Location’,’NorthEast') 

xlabel('k') 

ylabel('w(k)') 

axis([-pi()/a  pi()/a  0  3.5]) 


29 


Intentionally  left  blank 
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Appendix  C.  MATLAB  Code:  Two-Dimensional  Phonon  Dispersion  Curve 


function  Dispersion2D(epsilon,  sigma,  r,  m) 

%Phonon  Dispersion  Curves  for  Two-Dimensional  Monatomic  Lattice 
%function  Dispersion2D(epsilon,  sigma,  r,  m) 

%Calcualte  Interaction  Energy 

phi=  8  *  epsilon*  (((7 8  *  sigmaA  1 2)/ (rA  1 4))-((2 1  *  sigmaA6)/ (rA8))) . . . 

+8*epsilon*(((78*sigmaA12)/((2*r)A14))-((21*sigmaA6)/((2*r)A8))) 

%Next  nearest  neighbors  addition 

%Xenon  Test  Inputs 

%epsilon=0.02; 

%sigma=3.98; 

%r=4.33; 

%m=131.29 

syms  L  N  omega  k  kx  ky  a  b  K  G  M 

%Input  for  Monatomic  Homogeneous  Lattice: 

a=r; 

b=r; 

K=phi; 

G=phi; 

M=m*1.66053886e-27;  %convert  amu  to  kg 
%Lattice  Size: 

L=10;  N=10; 

%Set  loop  to  cover  entire  lattice 
e=N/2; 

j=i; 

for  i  =  -e:.l:e 
kx(j)=  (2*pi()/a)*(i/L); 
ky(j)=  (2*pi()/b)*(i/N); 

%k(i)=[kx(i),  ky(i)]; 

%Dispersion  Curve  Equation 

omega(j)  =  2*sqrt(K/M)*abs(sin((l/2)*kx(j)*a))... 

+2*sqrt(G/M)*abs(sin((l/2)*ky(j)*b)); 

j=j+i; 

end 

kx=double(kx); 

ky=double(ky); 

omega=double(omega); 

plot(kx,  omega) 

xlabel('k') 

ylabel('w(k)') 

axis([-pi()/a  pi()/a  0  1. 1  *max(omega)]) 
legend(’Acoustic  Branch', ’Location’, 'North') 
figure 

plot3(kx,  ky,  omega) 
xlabel('k');  ylabel('k');  zlabel('w(k)') 


31 


Intentionally  left  blank 


32 


Appendix  D.  Diamond  Structure  Control  File  (diamond.ctl) 


(set!  geometry-lattice  (make  lattice 

(basis-size  (sqrt  0.5)  (sqrt  0.5)  (sqrt  0.5)) 
(basisl  Oil) 

(basis2  101) 

(basis3  1  1  0) ) ) 

;  Corners  of  the  irreducible  Brillouin  zone  for  the  fee  lattice, 

;  in  a  canonical  order: 

(set!  k-points  (interpolate  4  (list 

(vector3  0  0.5  0.5)  ;  X 

(vector3  0  0.625  0.375)  ;  U 

(vector3  0  0.5  0)  ;  L 

(vector3  000)  ;  Gamma 

(vector3  0  0.5  0.5)  ;  X 

(vector3  0.25  0.75  0.5)  ;  W 

( vector3  0.375  0.75  0.375))))  ;  K 

;  define  a  couple  of  parameters  (which  we  can  set  from  the  command-line) 
(def ine-param  eps  11.56)  ;  the  dielectric  constant  of  the  spheres 
(def ine-param  r  0.25)  ;  the  radius  of  the  spheres 

(define  diel  (make  dielectric  (epsilon  eps))) 

;  A  diamond  lattice  has  two  "atoms"  per  unit  cell: 

(set!  geometry  (list  (make  sphere  (center  0.125  0.125  0.125)  (radius  r) 

(material  diel) ) 

(make  sphere  (center  -0.125  -0.125  -0.125)  (radius  r) 
(material  diel)))) 

;  (A  simple  fee  lattice  would  have  only  one  sphere/object  at  the  origin.) 

(set-param!  resolution  16)  ;  use  a  16x16x16  grid 

(set-param!  mesh-size  5) 

(set-param!  num-bands  5) 

;  run  calculation,  outputting  electric-field  energy  density  at  the  U  point: 
(run  (output-at-kpoint  (vector3  0  0.625  0.375)  output-dpwr) ) 
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Appendix  E.  Tri-Rods  Structure  Control  File  (tri-rods.ctl) 


(set!  num-bands  8) 

(set!  geometry-lattice  (make  lattice  (size  1  1  no-size) 

(basisl  (/  (sqrt  3)  2)  0.5) 

(basis2  (/  (sqrt  3)  2)  -0.5))) 

(set!  geometry  (list  (make  cylinder 

(center  000)  (radius  0.2)  (height  infinity) 
(material  (make  dielectric  (epsilon  12)))))) 

(set!  k-points  (list  (vector3  000)  ;  Gamma 

(vector3  00.50)  ;  M 

(vector3  (/  -3)  (/  3)  0)  ;  K 

(vector3  000)))  ;  Gamma 

(set!  k-points  (interpolate  4  k-points)) 

(set!  resolution  32) 


(run-tm  (output-at-kpoint 
(run-te) 


(vector3  (/  -3)  (/  3)  0) 

f ix-ef ield-phase  output-efield-z) ) 
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Appendix  F.  An  Introduction  to  Ferroelectric  Fatigue 


Introduction 

Ferroelectric  materials  have  many  applications  because  of  their  ability  to  have  spontaneous 
polarization.  They  are  used  in  sensors  and  actuators  in  addition  to  dynamic  and  non-volatile 
ferroelectric  memory  devices.  Unfortunately,  there  are  several  drawbacks  to  ferroelectric 
materials.  Ferroelectrics  have  low  fracture  toughness  and  can  crack  easily  because  of  their  brittle 
nature.  Another  main  issue  preventing  the  widespread  commercialization  of  ferroelectric  materials 
is  ferroelectric  (polarization)  fatigue.  This  is  the  loss  of  switchable  polarization  in  a  material  as  a 
function  of  the  number  of  switching  cycles.  A  switching  cycle  is  when  an  applied  external  electric 
field  is  cycled  from  one  direction  to  the  complete  opposite.  The  crack  fonnation  and  ferroelectric 
fatigue  lead  to  a  loss  in  performance  of  the  material.  Ferroelectric  performance  is  the  ability  and 
quickness  with  which  the  material  can  switch  polarization.  This  performance  loss  generates 
reliability  concerns  for  devices  that  use  ferroelectrics. 

A  way  to  predict  ferroelectric  fatigue  would  provide  great  insight  into  improving  the  fatigue  life. 
Fatigue  life  is  the  number  of  cycles  until  the  material  fails  and  can  no  longer  switch  polarizations. 
This  would  open  the  door  for  optimization  of  ferroelectric  materials.  One  way  to  accomplish  this 
goal  is  through  the  use  of  a  model  that  can  provide  fatigue  life  based  on  entered  parameters.  Arias 
et  al.6  propose  a  phenomenological  cohesive  model  of  ferroelectric  fatigue  to  predict  fatigue  life. 

Background 

A  dielectric  material  is  one  that  is  non-conducting.  Ferroelectrics  are  a  subset  of  dielectrics  that 
exhibit  interesting  properties  such  as  a  high  dielectric  constant.  They  are  materials  that  are  polar 
by  nature.  They  are  polar  because  ferroelectric  crystal  cells  have  a  dipole  moment  even  when  not 
subjected  to  an  external  electric  field.  The  dipole  is  a  result  of  different  center  locations  of  the 
positive  and  negative  charges  within  the  crystal  cell.  Because  of  the  dipole  moment,  the  crystal  is 
said  to  be  polarized,  meaning  it  contains  a  dipole  moment  that  acts  over  some  volume.  Crystals 
that  have  spontaneous  polarization  are  called  pyroelectric  and  the  polarization  occurs  along  the 
polar  axis.  Ferroelectrics  have  multiple  equilibrium  states  of  the  spontaneous  polarization  vector. 
An  applied  external  electric  field  can  switch  the  orientation  of  the  spontaneous  polarization  vector 
is  if  is  strong  enough  to  overcome  the  polarization  energy  barrier. 

Ferroelectric  materials  have  a  temperature  above  which  they  become  paraelectric,  meaning  they 
no  longer  are  spontaneously  polarized.  The  temperature  where  this  occurs  is  the  Curie  point,  Tc. 
Below  the  Curie  temperature,  the  ferroelectric  crystal  undergoes  a  structural  phase  change.  This 


6 Arias,  I.;  Serebrinsky,  S.;  Ortiz,  M.  A  Phenomenological  Cohesive  Model  of  Ferroelectric  Fatigue.  Acta  Mater. 
2006,  54  (4),  975-984. 
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shift  of  atoms  within  the  cell  generates  the  polarization  as  well  as  new  unit  cell  dimensions.  The 
spontaneous  polarization  does  not  have  to  occur  uniformly  across  the  whole  medium.  Multiple 
polarization  vectors  can  exist  across  the  material.  Areas  of  the  same  polarization  direction  are 
known  as  ferroelectric  domains.  The  orientation  of  ferroelectric  domains  throughout  the  material 
is  attributable  to  the  mechanical  and  electrical  boundary  conditions  to  which  the  material  is 
subjected.  The  domains  orient  themselves  to  minimize  the  surface  charge  and  the  elastic  energy 
resultant  of  mechanical  constraints.  The  surface  charge  is  a  result  of  the  onset  of  spontaneous 
polarization.  This  surface  charge  produces  an  electric  field  called  a  depolarizing  field  which  acts 
opposite  the  overall  polarization7.  The  area  where  two  domains  meet  is  a  domain  wall.  If  the 
domains  are  oppositely  polarized,  it  is  a  180-degree  domain  wall.  If  the  domains’  polarization 
vectors  are  perpendicular  to  one  another,  it  is  a  90-degree  domain  wall. 

A  plot  of  polarizations  versus  electric  field  will  produce  a  hysteresis  loop.  This  demonstrates  the 
ability  of  an  electric  field  to  reverse  the  spontaneous  polarization  of  some  or  all  domains  in  the 
material.  If  the  crystal  starts  with  an  equal  number  of  opposing  domains,  its  overall  polarization  is 
zero.  If  an  electric  field  is  applied,  some  of  the  domains  opposing  the  field  will  reverse.  As  the 
field  is  increased,  the  domains  will  continue  to  switch  until  a  critical  point  (the  state  of  saturation) 
is  reached.  At  this  point,  all  domains  are  oriented  in  the  direction  of  the  field.  If  the  field  is 
decreased,  the  polarization  will  not  return  to  zero  but  will  decrease  more  gradually.  When  the  field 
is  reduced  to  zero,  some  of  the  domains  will  retain  the  direction  of  the  first  critical  state;  this  is 
known  as  remnant  polarization,  Pr.  If  an  opposite  field  is  applied,  the  domains  will  reverse  direc¬ 
tion  until  all  are  oriented  with  the  field  at  a  second  opposite  state  of  saturation.  The  coercive  field, 
Ec,  is  the  electric  field  required  for  the  polarization  to  equal  zero. 

Ferroelectric  fatigue  is  a  very  complex  phenomenon.  It  is  defined  as  the  loss  of  remnant  polariza¬ 
tion  as  a  function  of  the  number  of  switching  cycles.  In  addition  to  electric  field  cycling,  other 
factors  can  decrease  the  remnant  polarization.  Several  fatigue  mechanisms  have  been  identified  as 
causes  of  ferroelectric  fatigue.  Domain  wall  pinning  is  where  space  charges  or  injected  charges 
prevent  the  movement  of  domain  walls;  this  mechanism  operates  mostly  in  the  bulk  of  the 
material.  Near  the  electrode  interfaces,  nucleation  of  oppositely  oriented  domains  is  inhibited. 

This  means  that  the  formation  of  opposite  domains  is  prevented  by  space  charges  produced  during 
electric  field  cycling.  Oxygen  vacancies  in  the  crystal  could  also  be  a  fatigue  mechanism.  The 
vacancies  could  pin  domain  walls  and/or  create  structural  damage  at  the  interface  between  the 
ferroelectric  and  electrode.  Micro-cracking  is  also  identified  as  a  cause  of  polarization  loss, 
although  it  is  a  purely  mechanical  defect.  Concentrated  stress  or  electric  fields  at  the  crack  tip 
induce  local  switching  that  can  shield  or  promote  crack  growth.  Crack  growth  is  affected  by 
electrical  conditions  at  crack  faces,  by  the  grain  size,  and  by  the  porosity  of  the  material8.  Many 


7Damjanovic,  D.  Ferroelectric,  Dielectric  and  Piezoelectric  Properties  of  Ferroelectric  Thin  Films  and  Ceramics. 
Rep.  Prog.  Phys.  1998,  61  (9),  1267-1324. 

o 

Arias,  I.;  Serebrinsky,  S.;  Ortiz,  M.  A  Phenomenological  Cohesive  Model  of  Ferroelectric  Fatigue.  Acta  Materialia 
2006,  54,  975-984. 
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fatigue  mechanisms  may  depend  on  frequency  as  well  as  the  amplitude  of  the  applied  electric  field. 
Ferroelectric  fatigue  causes  the  polarization  versus  electric  field  hysteretic  loop  to  change  shape. 
After  a  certain  number  of  cycles  corresponding  with  onset  of  fatigue,  the  maximum  polarization 
level  or  state  of  saturation  decreases.  The  crystal  is  no  longer  able  to  attain  its  original  maximum 
polarization. 

Cohesive  Model 

Paris  et  al.9  laid  the  foundation  for  fatigue  life  prediction.  The  work  provided  phenomenological 
laws  relating  applied  stress  amplitude  to  fatigue  crack  growth  rate.  The  model  is  known  as  Paris’ 
Law.  The  law  was  ground  breaking  at  the  time  but  was  based  on  the  ideal  conditions  of  small 
yielding,  constant  stress  amplitude  and  long  cracks.  If  these  criteria  were  not  met,  the  law’s 
prediction  capabilities  were  greatly  diminished.  Although  Paris’s  law  only  worked  for  select 
cases,  it  paved  the  way  for  future  fatigue  modeling  in  materials. 

A  model  using  cohesive  laws  was  developed  by  Nguyen  et  al.10  to  predict  fatigue  life  based  on 
fatigue  crack  growth.  The  model  was  prepared  to  predict  the  fatigue  life  of  metals  (aluminum  was 
used  for  experimental  data)).  The  model  augmented  Paris’  law  through  the  use  of  cohesive  laws. 
These  laws  exhibit  unloading-loading  hysteresis  which  allows  for  steady  crack  growth  modeling. 
Cohesive  theories  take  fracture  as  a  gradual  process  where  incipient  material  surface  separation  is 
resisted  by  cohesive  tractions.  The  tractions  go  to  zero  at  a  point  of  critical  opening  displacement. 
Previous  models  used  an  elastic  cycle  of  fatigue  but  the  cohesive  law  replaces  that  with  unloading¬ 
loading  hysteresis.  Many  of  the  previous  modifications  of  Paris’  law  had  been  ad  hoc,  meaning 
they  were  derived  from  experimental  data.  The  cohesive  laws  provided  a  means  to  create  a  model 
that  could  predict  fatigue  life  before  experimental  data  were  taken.  The  crack  growth  results  from 
interaction  of  bulk  cycle  plasticity,  closure,  and  gradual  decohesion  at  the  crack  tip.  This  model 
was  a  step  in  the  right  direction  but  no  cohesive  model  existed  for  ferroelectric  materials. 

In  2004,  Arias  et  al.  developed  a  model  to  describe  fatigue-crack  nucleation  and  growth  for 
ferroelectric  materials  subject  to  electro-mechanical  loading.  This  was  the  first  model  based 
around  a  hysteretic  cohesive  law  that  couples  mechanical  and  electric  fields  to  predict  fatigue  in 
ferroelectrics.  This  was  important  because  ferroelectric  materials  demonstrate  electrical  and 
mechanical  fatigue  under  cyclic  electrical  loading.  The  model  also  could  take  into  account  electro¬ 
mechanical  loading  which  is  what  ferroelectrics  are  subjected  to  in  operation  more  often  than 
purely  mechanical  stress  as  in  metals.  The  cohesive  law  can  be  used  in  tandem  with  general 
constitutive  relations  of  bulk  behavior  to  predict  fatigue  crack  growth  under  arbitrary  loading 
conditions.  This  model  is  also  able  to  predict  fatigue  crack  nucleation  which  the  previous  crack 
growth  model  (2001)  could  not.  This  is  important  for  fatigue  prediction  in  smooth-surfaced 


9Paris,  P.C.;  Gomez,  M.P.;  Anderson,  W.P.  A  Rational  Analytical  Theory  of  Fatigue.  The  Trend  in  Engineering 
1961,75,9-14. 

10Nguyen,  O.;  Repetto,  E.  A.;  Ortiz,  M.;  Radovitzky,  R.  A  Cohesive  Model  of  Fatigue  Crack  Growth.  Int.  J. 
Fract.  2001,  110  (4),  351-369. 
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material  with  no  initial  cracking.  One  large  benefit  of  this  model  is  the  ability  of  cohesive  theories 
to  apply  to  arbitrary  crack  and  specimen  geometries  and  loading  histories.  A  shortcoming  of  the 
model  is  that  it  is  phenomenological  in  nature  and  some  details  were  chosen  because  of  conveni¬ 
ence.  We  must  determine  some  parameters  of  the  model  experimentally  by  fitting  to  fatigue  data 
or  by  measuring  parameters  individually. 

The  most  recent  phenomenological  cohesive  model  of  ferroelectric  fatigue  produced  by  Arias  et  al. 
(2006)  is  a  great  advancement  in  ferroelectric  fatigue  life  prediction.  Like  the  preceding  model,  it 
predicts  ferroelectric  fatigue  under  electro-mechanical  loading  through  the  use  of  cohesive  laws. 
The  revised  cohesive  law  couples  mechanical  displacement  and  electric-potential  discontinuity  to 
mechanical  tractions  and  surface-charge  density.  The  model  applies  to  fatigue  when  it  is  localized 
in  planar  surfaces  that  are  treated  as  cohesive  surfaces.  The  surface  is  modeled  as  a  ferroelectric 
material  by  representing  both  mechanical  opening  displacement  and  electric-potential  discontinu¬ 
ity.  The  cohesive  law  is  extended  to  represent  these  and  mechanical  tractions  and  surface-charge 
density  as  a  work-conjugate  pair.  To  retrieve  qualitative  data,  a  “Smith-Ferrante  monotonic 
envelope”  and  an  exponential-decay  law  of  loading-unloading  hysteresis  were  used.  The  model 
delivers  the  following: 

-Existence  of  a  threshold  field  for  the  onset  of  fatigue, 

-Dependence  of  threshold  on  the  applied  field  frequency, 

-Dependence  of  fatigue  life  on  the  amplitude  of  the  field, 

-Dependence  of  the  coercive  field  on  the  size  of  the  component  (size  effect). 

To  gauge  the  validity  of  the  system,  a  simple  test  configuration  provided  experimental  data  to 
compare  to  the  model  output.  The  test  configuration  consisted  of  an  “infinite”  slab  of  PZT  (lead 
zirconate  titanate  Pb[ZrxTii_x]03)  acted  upon  by  an  oscillatory  voltage  differential  across  the  slab 
without  other  stresses.  The  results  were  not  conclusive  but  they  did  indicate  planar-like  regions 
affected  by  cycling  may  lead  to  the  observed  fatigue  in  tetragonal  PZT. 

Conclusion 

To  accurately  determine  the  interactions  between  fracture,  deformation,  and  defective  structures 
that  cause  ferroelectric  fatigue,  a  physics  based  multi-scale  model  is  needed.  This  model  would  be 
a  huge  step  in  predicting  the  fatigue  life  of  ferroelectric  devices.  Until  that  is  produced,  the  current 
phenomenological  model  should  be  improved  upon  and  validated  by  experiments  so  it  will  be 
useful  for  future  ferroelectric  design.  Several  actions  can  be  taken  to  improve  the  current  model. 
The  first  is  to  run  the  same  simple  test  configuration  that  was  done  with  PZT  on  different  ferro¬ 
electric  materials.  Having  more  data  on  multiple  materials  will  aid  in  the  following  proposed 
actions.  Extension  and  calibration  of  the  model  could  improve  data  agreement.  The  results  could 
also  optimize  fit  to  the  experimental  data.  Assumptions  made  in  the  derivation  of  the  model  could 
be  studied  further  to  confirm  their  validity.  Another  area  of  concern  is  the  fact  that  the  model  only 
applies  to  changes  in  properties  that  lead  to  fatigue  in  planar  regions.  A  more  comprehensive 
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model  working  on  arbitrary  sections  of  material  would  be  very  helpful  for  specific  material 
application  fatigue  prediction.  Full  finite  element  calculations,  more  precise  bulk  material 
constitutive  relations,  and  the  study  of  the  cohesive  law  aspects  of  the  monotonic  envelope  and 
loading-unloading  law  should  be  pursued. 

Although  results  were  inconclusive  for  the  cohesive  model  compared  to  the  experimental  test 
configuration,  it  is  crucial  to  put  work  into  the  model  to  improve  it.  To  have  a  model  that  could 
accurately  predict  the  ferroelectric  fatigue  effect  of  an  electro-mechanical  loading  would  greatly 
aid  future  work  with  ferroelectric  materials  such  as  PZT.  The  properties  of  these  ferroelectric 
materials  make  them  crucial  in  the  improvement  of  current  technologies  and  the  implementation 
of  emerging  technologies.  Fatigue  aspects  are  especially  important  for  the  Army  where 
equipment  is  constantly  subjected  to  harsh  and  dynamic  operating  conditions. 
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